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Quantitative phase imaging has become a topic of considerable interest in the microscopy community. We 
have recently described one such technique based on the use of a partitioned detection aperture, which 
can be operated in a single shot with an extended source [Opt. Lett. 37, 4062 (2012)]. We follow up on 
this work by providing a rigorous theory of our technique using paraxial wave optics, where we derive 
fully three-dimensional spread functions for both phase and intensity. Using these functions we discuss 
methods of phase reconstruction for in- and out-of-focus samples, insensitive to weak attenuations of 
light. Our approach provides a strategy for detection-limited lateral resolution with an extended depth of 
field, and is applicable to imaging smooth and rough samples. 
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1. INTRODUCTION 

Quantitative phase imaging reveals optical path variations in 
almost transparent (and also specularly reflecting) objects, and 
thus provides useful information for biological or industrial ap¬ 
plications. Different strategies exist to obtain phase information 
[1], most of which are interferometric and rely on monochro¬ 
matic illumination. In the case of quasi-monochromatic illumi¬ 
nation, optical phase is not well defined and interferometric 
techniques only measure changes in optical phase relative to a 
self-reference provided, for example, by spatial filtering [2-4], 
or shearing [5-7]. Alternatively, gradients in optical phase can 
be inferred non-interferometrically based purely on intensity 
imaging, using, for example, the transport of intensity equation 
[8-10], or by directly measuring changes in the direction of flux 
density, or local light tilt [11-17]. We have recently introduced 
a method to measure local light tilts that is based on the use of 
a partitioned detection aperture [18]. By associating local light 
tilts with phase gradients and integrating these over space, we 
thus obtain quantitative images of phase. The method is passive, 
single-shot, provides high spatial resolution, and works with 
quasi-broadband, partially coherent illumination. A reflection 
version of the method is described in [19], as is a version imple¬ 
mented for closed-loop adaptive optic wavefront sensing [20]. 

The purpose of this work is to present an in depth theo¬ 
retical analysis of our partitioned aperture wavefront imag¬ 


ing technique that goes significantly beyond the cursory and 
mostly experimental expositions we provided in our previous 
reports [18,19]. In particular, our previous reports were based 
on smooth phase approximations and on strict assumptions re¬ 
garding the numerical apertures of both the illumination and 
detection optics of our system, which ultimately limited spatial 
resolution. In this paper, we generalize to phases that are not 
necessarily smooth (though they must remain small), and we 
relax our assumptions regarding numerical apertures, enabling 
access to improved spatial resolution. Finally, we extend our 
analysis to three dimensions by considering samples that can be 
out of focus, and describing a refocusing strategy conceptually 
similar to that used in light field [21,22] or integral [23] imaging. 
The goal of this paper is to provide a rigorous theoretical under¬ 
pinning to partitioned aperture wavefront imaging, which has 
been long overdue. In the process, we also discuss strategies to 
expand on its capabilities. 

2. PHASE IMAGING WITH A PARTITIONED APERTURE 

A schematic of our method is shown in Fig. 1. A phase sample is 
modeled as a thin transparent layer of variable optical thickness. 
The sample is trans-illuminated by an axially symmetric beam 
of light of uniform intensity, which is characterized at each point 
by a distribution of light rays filling a well-defined square illu¬ 
mination numerical aperture. After traversing the sample, the 
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light is imaged by a device that comprises a lens of focal length 
f e , a partitioned lens assembly of focal length fy and a camera, 
arranged in a 3/ imaging configuration. For simplicity, we as¬ 
sume that the lens assembly is masked such that each individual 
lens in the assembly appears square. 


NAi f c f e f d 



illumination aperture 


Fig. 1 . (a) Partitioned detection aperture microscope in transmis¬ 
sion mode, (b) Image of the illumination aperture on the face of 
the partitioned assembly of four lenses: dashed square identifies 
the aperture in the absence of a phase sample, the shifted solid 
square demonstrates the effect of light refraction by the sample, 
(c) Four off-axis lens assembly: the locations of optical axes of 
the four lenses with respect to the optical axis of the system are 
identified by vectors b z , i = 1..4. A larger aperture (shading) 
ensures the detection apertures are square. 

The purpose of the quatrefoil lens assembly is to simultane¬ 
ously project multiple separated images onto the camera (here 
four), each acquired from a different oblique direction. Local 
tilts of light rays caused by the sample thus lead to local intensity 
differences in the projected images, as can be understood from 
Fig. 1. In the absence of a sample, the square illumination aper¬ 
ture forms a centered image on the surface of the lens assembly, 
and one obtains four images of uniform and equal intensities 
at the camera. When a phase sample is present, the intensity 
distribution between the images changes. Indeed, phase slopes 
at a given sample point cause the light rays through that point to 
tilt and the image of the illumination aperture from that point to 
shift (Fig. lb). Correspondingly, the four images of the sample 
point projected onto the camera exhibit unequal intensities. If 
the detection aperture is large enough there is no vignetting of 
the illumination aperture image and one can uniquely retrieve 
the tilts from the intensity distribution of the four images. As 
we discuss below, the retrieved light tilts provide a measure 
of sample-induced phase gradients, and hence, ultimately, of 
phase. 


A. Partitioned aperture imaging 

The four off-axis lenses form four images at the detector plane. 
To begin, we consider light propagating through one of the 
lenses shifted from the optical axis of the system by vector b. 
Throughout this work we confine ourselves to the paraxial wave 
approximation [24] and assume, for simplicity, unit magnifica¬ 
tion (i.e. f e = The field at the detector (camera) plane is thus 


given by 

Ed(t, t) = J d 2 toe -^ b 

xCSE d (i-i 0 )E 0 (i 0r t), (1) 

where A = 1/k is the wavelength, the coordinate r = r z — b 
accounts for the shift of the image with respect to the position 
b of the detection lens, and Fo( r O/f) is the illumination field 
at time t at the focal plane (in the absence of a sample). The 
detection numerical aperture of the system is defined by that of 
an individual off-axis detection lens fd, and the coherent spread 
function of each lens, when centered, is given by 

CSF d (r)=d) 2 / d V C 2 "^'A d (r' ) (2) 

where A d (r) is the square aperture function of an individual 
lens. 

The intensity at the detector plane is defined as I d (r) = 
(Ed(r / t)Ej(r,t)) / where the average (...) is taken over a time 
scale much longer than the coherence time of the light. The 
intensity is expressed as a weighted sum of the mutual intensity 
/o(ro c , r Qd) over the focal plane 

h(r) = Jd\ cd \ d e- t2 ^%(r 0 c^ 

xCSF d (r-r 0+ )CSF^(r —r 0 _), (3) 

where Jo(roc/ r od) — (Fo( r o+/t)£o( r O-/t)) characterizes spatial 
correlations of the illumination field between pairs of points 
in the focal plane, and tq± = tq c d= r od /2. The off-axis shift of 
the detection lens is accounted for by the phase factor with the 
spatial frequency proportional to the oblique detection tilt angle 
b/f e . Eq. (3) serves as the basis for what follows. 

B. Partially coherent illumination 

In our optical system a transparent sample is trans-illuminated 
by a uniform beam of quasi-monochromatic light resulting from 
Kohler illumination with a condenser lens of focal length f c . 
The spatial extent of the light source is defined by the aperture 
function A c (r). We consider a spatially incoherent source de¬ 
scribed by the function / c (r c ,r d ) = k~ 2 Io A c (r c )<5(r d ), where 
J 0 is the uniform intensity of the source, and J(r) is the two- 
dimensional Dirac delta function. The illumination mutual in¬ 
tensity at the focal plane (in the absence of a sample) is then 
given by Jo(r c , r d ) = Io^o( r rf)/ where the illumination coherence 
function is given by 

Mt) = /fk / d2r ' e ~ i27I ^ lA <^' < 4) 

and we defined the solid angle O c = f d 2 r|A c (r')| 2 //«? of 
the source as viewed through the condenser lens. As noted 
above, we assume a square illumination aperture of size 2 d c x 
2 d c . The illumination coherence function is thus separable 
^o(r) = sine (nx/af) sine (ny/af), where a\ = A/(2NAj) de¬ 
fines a characteristic length scale of the illumination mutual 
intensity, NAi = d c / f c is the numerical aperture of the illumi¬ 
nation, and sinc(x) = sin(x)/x. The solid angle of the source is 
Q c = (2NA;) 2 . 
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C. Thin phase sample 

A thin phase sample is characterized by a transmissivity function 
T(r) = e i( P( r \ The complex-valued phase of the sample (p(r) = 
cp(r) + ioc( r) accounts for spatial variations in the optical path 
and attenuation of light fields in the sample. The phase can be 
approximated by <p( r) = (2zr/A) J Q a dzn{t,z), where n(r,z) is 
the local index of refraction, and the thickness d( r) of the sample 
is assumed small. In the reflection mode, the real part of the 
phase is related to the topographic profile of the sample [19]. 
The field after propagation through the sample is E( r, f)T(r) = 
E( r, £)e l '^( r ). Correspondingly, the mutual intensity immediately 
after the sample depends on the relative phases between pairs 
of points 

Js{ r c ,r d ) = J(r C/ r d ) exp [z'</> (r c + y) - i<p* (r c - y)] , (5) 

where J(r c , r^) is the mutual intensity of the illumination field 
E( r, t) immediately before the sample. 


A. Imaging smooth phase gradients using a quatrefoil parti¬ 
tioned aperture 

The integration in Eq. (7) can be carried out explicitly when the 
illumination aperture A c (r) is a square. We assume that the 
coordinate system is centered at the optical axis of the system, 
and that the axes are aligned with the sides of the aperture. 
In a quatrefoil assembly, the centers of four lenses are placed 
symmetrically at positions 

bi = (b, b), b 2 = (-27, b) r b 3 = (-2?, -b), b 4 = {b, -b), (8) 

where b is the distance of a lens center from the two axes (x,y) 
(see Fig. lc). The four lenses of the assembly collect different 
amounts of light depending on the local tilt angles imparted on 
the light by the sample, and form images of varying intensity. 

Assuming the detection aperture of each lens exceeds the 
illumination aperture, from Eq. (7) we obtain 


[l± e * (r) l 

[uV r) l 

NAi _ 

NAi _ 


3. APPROXIMATION OF SMOOTH PHASE GRADIENTS 

The in-focus imaging of smooth phase gradients has been previ¬ 
ously described in Ref. [11, 16, 18] using geometrical optics. In 
this section we provide an alternative derivation of the formal¬ 
ism based on the wave optics, and identify its limitations. 

When we employ the smooth-phase approximation cp(r c ± 
r^/2) ~ (p(r c ) =b (r^/2) V(p(r c ) , the mutual intensity (5) at the 
focal plane immediately after the sample becomes 

Jo( r Oc/ r orf) « kM r od) exp [%V«p(r 0c ) - 2a(r 0c )]. (6) 

In this manner, sample-induced phase gradients become im¬ 
printed on the optical mutual intensity, provided the coherence 
function yo(r^) i s least partially coherent. 

The approximation (6) is valid provided a 2 \\7 3 (p\ <C |V<p| 
and aj\\7 2 oc\ <C oc, that is provided phases are smooth. When 
more stringent conditions +-|V 2 <p| <C |V<p| and + | Vx| <C oc 
are satisfied, that is when the phase gradients are smooth, 
we can approximate V<p( tq c ) « V<p(r c ) and oc(xq c ) « oc(r c ) 
in Eq. (6). Substituting the mutual intensity in Eq. (3), and 
also using the relation (4) between the Fourier components 
of the mutual intensity and the aperture function of the con¬ 
denser (see Appendix for our convention for Fourier transforms), 
fo(k) =A c (-/ c kA)/(fc 2 n c ) , we obtain the intensity at the de¬ 
tector plane 



This equation can be readily interpreted from geometrical optics. 
According to Eq. (7), light rays propagating within the cone 
limited by the illumination numerical aperture are tilted by a 
common angle at a given sample point. This tilt is translated 
into an effective lateral shift of the illumination aperture that is 
proportional to the tilt. The light is collected by the detection 
aperture A^, itself offset from the optical axis by vector b. Finally, 
the image is formed as an incoherent sum of the rays. The 
intensity is reduced exponentially according to the local value 
of attenuation parameter oc( r). 


where we have introduced the local tilt angles 0 X and 0 y along 
x and y axes, respectively. The connection between these local 
tilt angles of light and the local phase gradients of the sample is 
given by 

Ox = 2 —d X <P/ Oy = ——d y (p, (10) 

which constitutes the general crux of intensity-based phase imag¬ 
ing techniques [25]. We note that the signs of tilt angles in Eq. (9) 
are defined by the quadrant of the detection lens. For example, if 
the lens occupies the fourth quadrant as seen from the incoming 
beam, that is b = (+27,-2?), one should take the signs (+) and 
(—) for the x and y components of the tilts, respectively. The fac¬ 
tor 1/4 in Eq. (9) accounts for the splitting of the lens assembly 
into four equal parts. 

The light tilt angles are extracted from the four images at the 
camera using the linear combinations: 


6 X = NAi (h + h-I 2 -h)/I to t, 

0 y = NA { (h + h-h-k)/hot, (ID 

where 4 for k = 1..4 are the intensities detected in the four 
quadrants of the detector plane, and the total intensity is If 0 t = 
Ek=l.A 4- According to Eq. (9), the range for light tilt measure¬ 
ments is defined by the illumination numerical aperture, such 
that 1 0 X/ y | < NAi. The tilt angles can be expressed using the 
wavevector components, 0 x>y = k X/ y/k , so that the range of the 
tilts can be also written as \k x ,y\ < l/(2+). Eqs. (10) and (11) 
have been previously derived in Ref. [11, 16, 18] using geometri¬ 
cal optics. 

Eqs. (10) and (11) can be used to reconstruct the sample 
phase distribution. Assuming that the spatial support of 
the sample phase cp( r) is finite within the imaging field of 
view, we find the Fourier components of the tilts (11) and 
use Eq. (10) to obtain the spectral representation of phase [6] 
<p(k) =# —ik [0 X ( k) + i0y( k)] /(k x + iky). The phase profile is 
found by using the inverse Fourier transformation of <p(k) 


*/ 


(r) = -ik / d z ke l 


j2i, Jlnkr ^4(k) + i0 y ( k) 


k x + iky 


( 12 ) 


where the integration extends over the bandwidth of optical 
system limited by the detection aperture. The reconstruction (12) 
was derived in [18] without mention of the spatial resolution 
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of the method. Since the derivation assumed that phase gra¬ 
dients are smooth on the spatial scale of the mutual intensity, 
we conclude that (12) is accurate for spatial frequencies of the 
phase much smaller than the tilt range, \k x ,y\ -Cl/ (2 aj). In 
the following section we derive an alternative phase reconstruc¬ 
tion method characterized instead by a detection-limited spatial 
resolution. 


4. APPROXIMATION OF SMALL PHASES 

We have derived intensity (7) within the approximation of 
smooth phase gradients. We now generalize our formalism 
to the approximation of small phases, and allow for the pos¬ 
sibility of non-smooth phase distributions. In the process, we 
also extend our formalism to three dimensions by allowing the 
possibility of sample defocus. 

We assume that a phase object is located some distance z away 
from the focal plane (z > 0 correspond to displacements away 
from the imaging lens f e ). Within the approximation of small 
phases, the mutual intensity (5) at the sample plane becomes 

Js(ic,id) « J(ic,id) [1 + i<P (r+) - i<P* (r-)], (13) 

where r± = r c dt r d /2, and /( r C/ r d ) is the mutual intensity im¬ 
mediately before the sample. The small expansion parameter 
is the phase difference between two points within the spatial 
range of the mutual intensity. 

We propagate the illumination mutual intensity from the fo¬ 
cal plane toward the sample by the off-focus distance z, where it 
is modulated according to Eq. (13). The function is then propa¬ 
gated back to the focal plane by distance — z, where it reads 

Hrodod) « I 0 (^) 2 /d 2 r c d 2 r dW (r rf ) e -^‘^-^)^-^) 

X [1 + i(p (r+) - icp* (r_)], (14) 


where r± = r c ± r d /2. 

Substituting (14) into Eq. (3), we find that the intensity at the 
detector plane contains two terms 


Jd(r,z) = 4(r,z) + I 9 { r,z). 


(15) 


The first term I a accounting for weak attenuation of light in 
the sample is given by 


4( r,z) = I 0 J d 2 t 0 P x (t-i 0 ,z) 


\ -«( r o) 


(16) 


Using the notation q + = q + k/2 and q_ = q — k/2, the 
intensity transfer function reads 


< (k ' 2) = pk/ d 


2^ gQTTjkq 


fe 


A d fq + -b A* fq__ b 


fe 


fc 


Ac — U + +A c M-^q_ 


fc 


and the phase transfer function reads 


PJk,z) = - 


k 2 fX 


/ 


d 2 q e‘ 


,z 27 r|kq 


xA d f q+ _ b AUfq_-b 


Ac -^q + "A* -£q_ 


(19) 


( 20 ) 


These functions are calculated in Appendix in the case when 
both illumination and detection apertures are square in shape, 
and the numerical aperture of illumination does not exceed 
that of detection, NAi < NA^. Similar analysis of the transfer 
functions describing other phase-sensitive optical systems has 
been presented in Refs. [26-28]. 

Illustrations of the spread functions of phase and intensity 
are shown in Fig. 2 for the specific case NA^ = 2NAi, and for the 
in-focus position z = 0 (see Appendix for details). We observe 
that the spatial resolution associated with the transfer functions 
is a d = A / (2NAd), determined by the detection numerical aper¬ 
ture NAd (also see below a detailed discussion of the transfer 
functions). The phase spread function vanishes when averaged 
over the detection plane, indicating that it is sensitive to phase 
variations only. The average of the intensity spread function is a 
constant f (i 2 rP a (r,z) = 1/2 independent of axial position z. 


(b) 

* 

(a) 

(f) 

* 

(e) 

X 

(c) 

* 

(d) 

* 

(g) 

(h) 

* 


--i i 

-2 0 2 Position, units of a d 


where we refer to P a (r, z) as the intensity spread function. When 
attenuation of light is absent (oc = 0), the intensity becomes a 
constant independent of out-of-focus position z. 

The phase term is given by 

I<p( r,z) = I 0 J d 2 r 0 P,p(r-r 0 ,z)f(ro), (17) 


Fig. 2. (a)-(d) Spatial dependence of the phase spread functions 
a^P(p(r,z = 0) in quadrants 1-4, respectively; (e)-(h) the intensity 
spread functions fl 2 P a (r,z as 0) in quadrants 1-4, respectively. 
The functions are calculated for the in-focus position z = 0 and 
detection aperture NA^ = 2NAi. Lateral positions are in units 
of a d . 


where P^ (r, z) is the phase spread function. 

Concise representations of the intensity and phase spread 
functions are obtained using Fourier transforms. Assuming 
finite spatial support of <j>( r), we find the spatial frequency com¬ 
ponents of intensity (15) 

4(k)/Io = [S( k)/2 - a(k)] P a (k,z) + P^(k,z)^(k), (18) 

where P a (k,z) and P^(k,z) are the intensity and phase transfer 
functions respectively. 


Eqs. (19) and (20) are some of the main results of this work. 
Though we have taken as examples square apertures here, they 
are valid for arbitrary shapes of the detection and illumination 
apertures and for arbitrary defocus positions z. 

A. Phase-gradient expansion 

In the limit of small wavevectors k-^0we recover the approxi¬ 
mation of smooth phase gradients (7). Indeed, since our optical 
system can only detect variations of phase, the zero frequency 
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term of Eq. (20), corresponding to a constant phase, vanishes, 
P(p (k = 0,z) = 0. The first non-vanishing term is obtained 
by Taylor expanding (19) and (20) at small spatial frequencies, 
leading to 


h(?) 



X 1 


2*M 


fc 

Ink 


V r ?(r) 



A c (-r') , (21) 


which, in fact, is the Taylor expansion of Eq. (7) in the limit 
of small phase gradients (see details in Appendix). Specifically, 
when a square illumination aperture is smaller than the detection 
aperture of each lens of a quatrefoil assembly (8), we obtain the 
small phase-gradient approximation of Eq. (9) 


h M « j [i - 2a(r) ± ± ^3 y <p(r)] . (22) 


The phase gradients are orthogonal here to the edges of illumina¬ 
tion aperture. As before, the signs of phase gradients are defined 
by the detection lens quadrant. Phase gradients are obtained 
using the linear combinations of intensities (11). 


5. IN-FOCUS PHASE IMAGING 

The range of spatial frequencies where Eq. (12) remains valid 
can be estimated by studying the structure of the phase transfer 
function (20) beyond the Taylor expansion. We assume that a 
phase sample is located at the focal plane z = 0 and consider a 
square illumination aperture that is smaller than the detection 
aperture of each lens in a quatrefoil assembly (8). As an example, 
we assume that the detection aperture is a square of size 2 dp x 
2 dp. The corresponding numerical aperture is given by NA d = 

fe- 

Having derived the transfer functions (19) and (20), we can 
now refine our phase imaging method (12). Previously we made 
use of Eq. (10) to relate local light tilts to sample phase variations. 
Here, we formally define light tilts according to Eq. (11) and 
derive the relations between the tilts and phase 


0*( r ) 

= NAi J d 2 r 0 P e *(r-ro)<p(ro), 


0y(r) 

= NAj J d 2 r 0 Pg (r — r 0 )<p(r 0 ). 

(23) 


The spread functions of the tilts are linear combinations of the 

spread functions of intensities P^ (r), where m = 1..4 identifies 
the quadrant of the lens assembly: 

px _ n( 1) . p( 4 ) _ p(2) _ p(3) 

1 9 — 1 cp i 1 cp 1 cp 1 (p / 

P y e = P^+P^-P^-P^. (24) 

The transfer functions of the tilts Pg (k) and Pg (k) at the in¬ 
focus position z = 0 are obtained by substituting Eq. (15) in (11) 
and using the Fourier transform: 

&(k) = NAjPg (k)<p(k), 0 y (k) = NAjPg (k)<p(k), (25) 

We substitute Eq. (20) into Fourier-transformed Eqs. (11) and 
derive separable representations 

PeW =2ig-(k x )g+(ky), P^(k) = 2ig-(k v )g+(k x ), (26) 

where the functions g-(q) and g+(q) are illustrated in Fig. 3 and 
explicitly expressed in Appendix. 



Fig. 3. Plots of g- and for in-focus imaging when the detec¬ 
tion aperture is a square and NA d = 2NAi. 


We can understand the behavior of these functions qualita¬ 
tively, again from Fig. 1(b). Specifically, the function g-(q) is 
proportional to the tilt angle along one direction. Let us consider 
a prism-like sample that refracts light along y- axis. A positive 
tilt is imposed by the prism 6 y = k y /k > 0 along y- axis. When 
the tilt is relatively small, the image shifts upwards, and the two 
upper lenses of the quatrefoil assembly collect more light than 
the lower ones. Correspondingly, function ( k y ) grows linearly 
with k y until the tilt angle reaches the limiting value of the tilt 
range, 6 y = NA{. For larger angles, NAi <0y <NA d -NAi, 
the lower lenses collect no light, while the upper ones collect the 
maximum power. The power does not change until the image of 
the illumination aperture reaches the top edge of the lens assem¬ 
bly at 6 y = NA d — NAi. Still at larger angles of refraction, the 
intensity collected by the upper lenses begins to decrease until 
it falls to zero when the image of illumination aperture leaves 
the lens assembly altogether at 6 y = NA d . Similar arguments 
apply to negative tilts 6 y < 0, and it is straightforward to see that 
g- (k y ) is an odd function of k y . The optical bandwidth of the 
system is defined here by the detection aperture, | k x , y I < 1/flrf/ 
where a d = A/(2NA d ). The function g_ is shown in Fig. 3. 



Fig. 4. Plots of g- and for in-focus imaging when the detec¬ 
tion aperture is a square and NA d = NAi. 

The detection and illumination apertures completely charac¬ 
terize the light-tilt transfer functions. As we discuss below, to 
achieve detection-limited resolution for imaging out-of-focus 
phase objects, the numerical apertures of illumination and de¬ 
tection should be matched, NAi = NA d . Fig. 4 illustrates the 












































Research Article 


Journal of the Optical Society of America A 6 


transfer functions in this case. In contrast to Fig. 3, the range of 
frequencies for which =1/2 has been reduced to a point. 



Fig. 5. Spatial dependence of the tilt spread functions: 
fl?P^(r,z = 0) (a), and a?P^(r,z = 0) (b), calculated for the 
in-focus position z = 0, when the detection aperture is a square 
and NA d = 2NAi. Lateral positions are in units of a#. 


A. Spatial frequencies limited by the illumination aperture 

To begin, we consider spatial frequencies within the tilt range 
of our system, namely \k X/ y\ < 1/(2fl z ), as defined by the illumi¬ 
nation numerical aperture. The functions g- ( k Xf y ) = k X/ yUj and 
g+(k X/ y) = 1 — \k X/ ydj\ are then independent of the detection 
aperture. The corresponding Fourier components of the tilts (26) 
are given by 


S*(k) = 

iy( 1 -‘l Vzl)<P( k )' 


fly(k) = 

iy( 1 - M;|)<p(k), 

(27) 


functions to deconvolve light tilts and obtain the phase 




1 f c j2j ce ;27Tkr ftr(k) + ^y(k) 
NAii P*(k)+ijf(k)' 


(30) 


where integration now extends over the detection bandwidth. 
For a square detection aperture, the spread functions of light 
tilts are given in Appendix. Corresponding illustrations of the 
tilt spread functions P£ (r) and Pjj (r) are shown in Fig. 5 for 
the in-focus position z = 0. As can be observed, Pg (r) is even 
with respect to y -axis and odd with respect to x-axis, and for the 
function Pjj (r) vice-versa. 



-0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 

Position, units of a t 


Fig. 6. Reconstruction of a periodic phase profile cp(x)/cpo in 
Eq. (31) when Q = 1/aj and NA d = 2NAi using different meth¬ 
ods: Eq. (12) (solid blue line), Eq. (29) (red dashed line), and 
Eq. (30) (green solid line). The exact profile is depicted by the 
green solid line. Lateral position is in units of « z . 


where \k X/ y\ <1/ (2 a{). These equations uniquely map the direc¬ 
tions of light rays k X/ y /k to the tilt angles along the same axes 
in proportion to the Fourier components of the phase. These 
may be compared with Eqs. (10) utilized for our smooth phase 
approximation, and are found to be equivalent except for an 
additional attenuation of spatial frequencies orthogonal to the 
tilt directions. 

Phase is extracted from (27) using a linear combination of the 
tilts 


q>(k) = —ik 


8 X (k)+i8y(k) 

M1 - \kyUi\) + ik y {\ - \k x Oi\y 


(28) 


where \k x ,y\ < 1/(2a,), or, correspondingly. 


f(r) = -ik J 


d 2 ke J 


iinkr 


g*(k)+jgy(k) 


Ml - \ky*i\) + ik y (l ~ IMzl)' 


(29) 


where integration extends over the tilt range \k Xf y\ < 1/(2 U{). 
In this derivation we assumed that the spatial support of cp( r) 
is finite within the imaging field of view. We note that Eq. (29) 
is applicable to arbitrary shapes of the detection lens provided 
the associated aperture exceeds the numerical aperture of illumi¬ 
nation. The phase reconstruction (29) should be used when the 
illumination aperture is well characterized, while the detection 
aperture is not precisely known. However, the method does not 
account for the spatial frequencies beyond the tilt range. 


B. Spatial frequencies limited by the detection aperture 

The resolution of phase reconstruction can be improved if the 
detection aperture is precisely known. We use the known spread 


Let us compare the reconstruction methods (12), (29) and (30) 
in the simple case of a one-dimensional periodic phase profile 

cp(x) = (po cos(2ttQx), (31) 

where Q > 0 is the periodic spatial frequency. It is easy to verify 
that when Q < 1/ (2 af) different phase reconstruction methods 
yield an identical profile <p r ec(x) = (p(x). However, for larger 
frequencies the reconstructions are different. We consider for ex¬ 
ample frequencies within the range 1/(2 ap) < Q < l/a^ — l/a z -. 
The tilt transfer function is then constant (d=Q) = ±1/2 (see 
Fig. 3). The reconstruction results are summarized in Fig. 6, 
where NA d = 2NAi and frequency Q = 1 / a j. According to 
Eq. (12) we obtain <p r ec(x ) = ^cos(2ttQx), where the recon¬ 
structed amplitude ip = (po/ (2 UjQ) incorrectly depends on the 
spatial frequency Q instead of being the constant (po (Profile (I) 
in Fig. 6). This inconsistency is due to the incorrect mapping of 
phase to tilts outside of the tilt range. A refined approach (29) 
is only sensitive to spatial frequencies within the tilt range, so 
when Q > 1/ (2aj), the reconstructed profile (p r ec(x) = 0 (Pro¬ 
file (II) in Fig. 6) is also incorrect. Finally, Eq. (30) properly ac¬ 
counts for all spatial frequencies within the detection bandwidth 
Q < 1 / a dr an d, as a result, provides the correct reconstruction of 
the original profile (p r ec(x ) — f(x) (Profile (III) in Fig. 6). Thus 
all the spatial frequency components of an arbitrary phase pro¬ 
file within the detection bandwidth can, in principle, be found 
using Eq. (30). Approximate reconstructions of phase samples 
characterized by even broader spatial frequency distributions 
are also possible, as we discuss below. 
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The phase reconstruction methods (29) and (30) are the main 
results of this work. Experiments reported in Refs. [18,19] em¬ 
ployed the simplified method (12) of phase reconstruction. Since 
the spatial frequencies of the imaged samples were safely con¬ 
fined to the tilt range, the reconstruction (12) provided good 
approximations of phase profiles. 

C. Imaging a phase step 

It is instructive to apply the above formalism to the imaging of a 
phase, which manifestly is not smooth and features high spatial 
frequencies. The strong light diffraction from the discontinu¬ 
ity allows to visualize spatial resolution of our phase imaging 
method. For simplicity we consider a one-dimensional step 

f(i) = tp 0 H(x), (32) 

where the amplitude is small, |<pol •C 1, and the Heaviside 
function H(x) = 1 at x > 0, and H(x) = 0 at x < 0. 



Fig. 7. Detected and approximate light tilts 6 X / 0q as a function 
of position when imaging a phase step (32)and in the case of 
matched numerical apertures NA^ = NA{. Lateral position is in 
units of Uj. 

The Fourier components <p(k) are calculated using standard 
regularization, in which we add an infinitesimal imaginary part 
to the wave vector k x —> k x — iO 

« k > - ,33 > 

The Fourier components of the detected light tilts (11) resulting 
from the object are limited by the detection aperture 

&(k) = e 0 ^j^-S(k y ), 0y (k) = 0, \k x \ < 1 /a d , (34) 

Kx 

where we define the light tilt 6q = (cpo/ n)NA{. The disconti¬ 
nuity (33) diffracts light in all directions along x-axis, and the 
detection aperture collects only part of this light. 

In real space the detected tilt angle becomes 

/■l/flrf -r* 7 

0 x (x)=e o d k x e l2nkxX g-{k x )/k x , 0 y (r) = 0. (35) 

J-l/a d 

The dependence of the detected light tilt on the position along 
x-axis is plotted as a solid blue line in Fig. 7, assuming a square 
detection aperture of numerical aperture NA^ as NAi. For 
comparison we also plot an approximate tilt (dashed red line) 
given by 6 x (x) = 6q sinc(7rx/A z -), which is obtained from the 
detected tilt (34) by setting to zero the spatial frequencies | k x | > 



Fig. 8. Step profiles cp(x)/cpQ reconstructed using Eq. (12) (Profile 
I, red dashed line), and Eq. (29) (Profile II, solid blue line), with 
detection aperture NA^ = NA^. The profiles approximate a step 
function (32) shown as in solid black. Lateral position is in units 
of a d . 

1 / (2a,). We interpret the peak of the approximate tilt as follows. 
As a light ray refracts from the step, it effectively traverses an 
optical path h = cp^A/ (2 tc) within a lateral spot of size estimated 
by the resolution scale aj = A/(2NAi). The effective angle of 
refraction (i.e. tilt) is then estimated as the ratio of the optical 
path to the spot size, 6 = h/dj, which coincides with the peak 
value 6q = (cpo/n)NA{ of the approximate tilt. The detected 
tilt (shown in blue solid line) is obtained similarly using the 
detection-limited resolution. Therefore, the peak detected tilt is 
expected to be proportionally larger. However, since the transfer 
function g- in Fig. (4) is not linear, the peak value of the detected 
tilt is only somewhat larger than the peak of the approximate 
tilt. 

Results of phase reconstruction are summarized in Fig. 8. 
All reconstruction methods define the phase profiles up to an 
arbitrary constant, which we choose to be cp( 0) = cp o /2. The low- 
frequency part of the tilt 0 X ( k) = 0QajS(ky) valid within the tilt 
range \k x \ < 1/(2a z ), is employed in the reconstruction (29). The 
corresponding profile (I) is cp(x)/cpo = 1/2 + (l/7r)Si(7rx/fl z -), 
where Si(z) = J Q Z dq sin(q) / q, shown as a dashed red line. The 
oscillations near the discontinuity x = 0, known as the Gibbs 
phenomenon, are attributed to the finite range \k x \ < 1/ (2ai) 
of the integration in (29). The reconstruction (30) is given by 
cp(x)/cpo = 1/2 + (l/7r)Si(27rx/fl z ), and plotted as the solid 
blue line Profile (II). We note that in the considered case = a#. 
The oscillations around the step identify the spatial resolution 
of our system since they are a result of the finite detection band¬ 
width. Naturally, as the detection bandwidth is increased, the 
profile approaches the step function. 

6. OUT-OF-FOCUS PHASE IMAGING 

In this section we analyze phase imaging of out-of-focus objects 
located a distance z from the focal plane. The acquired images 
in this case exhibit both blurring and quadrant-dependent lat¬ 
eral shearing compared to in-focus imaging. As in light field 
imaging [22], we make use of the fact that shearing provides 
information on the defocus position and can be compensated 
for. The blurring of images for large defocus, however, cannot 
be compensated and degrades spatial resolution. We show that 
out-of-focus phase profiles can nevertheless, in principle, be re¬ 
constructed, though with a reduced spatial resolution depending 
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on the extent of defocus. 

The exact relation between the numerical apertures of illumi¬ 
nation and detection is not critical for phase reconstruction of 
in-focus objects, provided the detection aperture is accurately 
known. For defocused objects the situation is different, since the 
two apertures both control the effects of blurring and shearing. 
The case of matched numerical apertures NAi = NA d is special, 
since in that case lateral shearing observed in different quad¬ 
rants at the optical resolution uniquely identifies the defocus 
distance (see below and also Appendix). In the unbalanced case 
NAi < NA d , the situation is more complicated. The problem 
comes from the spectral components of intensities outside the 
tilt bandwidth, since they shear differently than the spectral com¬ 
ponents within the bandwidth. That is why, we consider only 
the frequencies within the detection bandwidth \k Xf y\ < l/(2a z ), 
which can be accomplished by applying a low-pass filter to the 
recorded images. The formalism presented below concentrates 
on the detection-limited phase reconstruction for a system with 
matched numerical apertures. However, it can also be applied 
in the unbalanced case after numerically filtering the bandwidth 
of the recorded images to be within the tilt bandwidth (see also 
Section H in Appendix). 



Fig. 9. Averaged phase spread functions, a z (P,y (x,y))y (solid 

( 2 ) 

blue line) and the negative of a z -(P^ ; (x,y))y (dashed red line), 
calculated in quadrants (1) and (2), as functions of lateral po¬ 
sition x at the defocus position z = z Z/ in the case of matched 
numerical apertures NA Z = NA d . The two functions are replicas 
of one another, shifted by Ax = 2NAiZ. Lateral position is in 
units of cij. 


The principle of numerical refocusing is illustrated in Fig. 9, 
where we considered the case of matched numerical apertures 
NAi = NA d . We find that defocused spread functions of quad¬ 
rants (1) and (2) averaged over the y coordinate, (P^ (x,y,z))y 

and (Pg, ’ (x, y, z)) y , are replicas of one another shifted along x- 
axis, up to a negative sign. Since the intensity of an image is 
obtained as a convolution (17) of phase and the phase spread 
function, we observe that the images of an arbitrary phase object 
shift accordingly. The relative shift is proportional to the defocus 
distance. Ax = 2NA|Z, and, thus, unambiguously identifies this 
distance. 

In the two-dimensional case, to avoid ambiguities, we also 
need to consider the first and fourth quadrants, and average 
the corresponding phase-sensitive intensities over the x coordi¬ 
nate. The averaged intensities depending on coordinate y are 
the shifted replicas of one another, up to a negative sign. The 


two relations between the intensities are given by 

(I<p^ ( x — NAiZ,y))y = —(ljp\x + NA\z,y))y, 

{I<p (x,y — NAjz))* = -<4 4 W + NAjz))*, (36) 

where (..) X/ y is the average over coordinates x and y, respectively, 
and I^ (r) is a phase-sensitive part of the intensity recorded in 

quadrant m. The intensity I^ (r) is obtained by subtracting the 
background from the total intensity (15). Identities (36) provide, 
in principle, an estimate of the defocus position z, for example, 
by maximizing the cross-correlation function of the averaged 
intensities. 

Light tilts adjusted for the shearing are given by 


e x (t,z) 


Oy{t,z) 


NAi h(x-,y-) + h(x-,y+) 

Itot -I 2 (x+,y-) - I 3 (x + ,y+) 


NAi h(x-,y-) + I 2 (x-,y+) 

hot -I 3 (x+,y-) - I i (x+,y+) 


(37) 


where x± = x d= NAiZ, y± = y db NAiZ, I m (r) is the intensity 
recorded in quadrant m, and If 0 t = Ii(x_,y_) + l 2 (x_,y+) + 
h{x+,y~) + h(x+,y+). 



Fig. 10. Intensities 1^ (x,z) and negative I® (x, z) as functions 
of lateral position x, calculated for different defocus distances 
and with matched numerical apertures NAi = NA d . Solid blue 
line and dashed red line correspond to z = z z ; solid green line 
and dashed magenta line correspond to z = 5z z . Lateral position 
is in units of a z . 


We find the Fourier components of (37) and use the known 
transfer functions of tilts (69) (see Appendix) to obtain the phase 
of an out-of-focus object 


<p(r,z) 


1 f gi2nkr ^x(K_z) + ffi/(k.z) 
NAj J P*(k,z) + zlf(k,z)' 


(38) 


where the integration extends over the detection bandwidth 
\k Xf y \ < 1/flj, when the numerical apertures of illumination and 
detection match, NA Z = NA d (i.e. a d = a z ). We note that for 
relatively small defocus, deconvolution is in principle possible 
within the detection bandwidth, based on the transfer func¬ 
tions (69). In cases when NAi < NA d the integration extends 
over the tilt bandwidth \k Xf y\ <1/(2a z ), which leads to a re¬ 
duced spatial resolution compared to the case of matched aper¬ 
tures. For a relatively large defocus, the combination of transfer 
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functions in (38) has zeros at k X/ y ^ 0. The bandwidth should be 
reduced to exclude these singularities from the integration. The 
appearance of singularities identifies the defocus regime where 
the blurring cannot be accurately compensated, which results 
in reduced resolution of phase reconstruction depending on the 
defocus distance. 

A. Imaging a defocused phase step 

We illustrate reconstruction (38) of a defocused phase step (32) 
by a system with matched numerical apertures NAi = NA^. 
Intensities detected in the first and second quadrants are ob¬ 
tained using the transfer functions (62). The characteristic axial 
length scale of defocusing is the depth of focus z z = A / (2NA 2 ). 
Intensities calculated at z = z z and z = 5z z are shown in Fig. 10. 
From these plots we can unambiguously identify the defocus 
positions. Indeed, pairs of intensities calculated at the axial 
positions z = z z and z = 5z z are shifted with respect to one 
another by distances bx = 2+• and bx = 10a z , respectively. Phase 
profiles reconstructed according to Eq. (38) are demonstrated in 
Fig. 11. We choose an arbitrary constant so that <p(0) = cpo/2 . 
The reconstructed profile cp(x) /cpo = 1/2+ (1 /n)Si(2nx/af) 
calculated at z = z z is identical to the in-focus one. One can 
show that detection-limited phase reconstruction is possible 
for relatively small defocus |z| < 2z z . In contrast, the profile 
reconstructed at z = 5z z is characterized by a reduced spa¬ 
tial resolution. The reason is that at this relatively large de¬ 
focus the transfer function Pq(Jc X/ z) = 0 at \k x \ — k c — 0.32/ +• 
within the detection bandwidth \k x \ < 1/+-. To avoid the sin¬ 
gularities in (38), we thus restrict the bandwidth of integration 
to a smaller range \k x \ < k c . The resulting profile is given 
by <p(x)/(po = 1/2 + (1 / n)Si(2nxk c ). Oscillations observed 
around discontinuity at x = 0 are due to the limited bandwidth 
of integration in (38). As the defocus distance is increased, the 
integration bandwidth is gradually reduced. We find that even 
at a significant defocus distance z = 20z z , the first zero of the tilt 
transfer function occurs at k c = 0.16/+, so that the bandwidth 
is reduced only two times compared to z = 5z z . 



Fig. 11. Phase profiles cp(x)/ cpo of a step (32), depicted by in 
solid black, reconstructed according to Eq. (38): the solid blue 
line corresponds to z = z z , dashed red line corresponds to z = 
5z z . Lateral position is in units of +•. 


7. CONCLUSION 

In this work we considered quantitative phase imaging us¬ 
ing a partitioned detection aperture microscope introduced in 
Ref. [18]. Our single-exposure method is characterized by a 


simple optical design that is light efficient and insensitive to po¬ 
larization. The microscope may be implemented in transmission 
or reflection modes [18,19]. Within the wave-optic paraxial ap¬ 
proximation we derived the three-dimensional spread functions 
of intensity and phase (19) and (20), analyzed their properties, 
and refined the phase-imaging approach of Ref. [18,19], making 
it suitable for imaging smooth and rough samples. The formal¬ 
ism is summarized in Eq. (30). In particular, we demonstrated 
that the spatial resolution of our method is limited by the detec¬ 
tion aperture in contrast to Refs. [18,19] where the resolution 
was defined by a smaller tilt range. We also showed how our 
method reconstructs the profile of a phase step provided is op¬ 
tical path-length is smaller than the wavelength. Further, we 
extended the method to imaging of defocused phase objects, as 
described in Eqs. (36), (37) and (38). Even for relatively large 
object defocus, the method provides precise reconstructions of 
rough and smooth samples with spatial resolution depending 
on the defocus distance. 
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9. APPENDIX 

We provide calculations of spread functions, which, for the sake 
of clarity, were omitted in the main text. 

In our calculations we employ the following convention for 
Fourier transforms 

/(r) = J d 2 k e i27rkr /(k), /(k) = J d 2 ie- i2nkl f{i). (39) 


A. Phase-gradient expansion of transfer functions 

In the limit of small wavevectors k —»• 0, one can simplify (20). 
Assuming that the aperture is a real function. A* = A c , we 
expand the difference of the functions into the Taylor series and 
keep the first non-vanishing term 

Ac q +) - A* (-f q-) «kV q A c (-=|q) , (40) 


where q± == q ± k/2. Substituting this expression in Eq. (20), 
and restricting the expansions to terms linear in k, we derive an 
approximate phase transfer function 


P<p( k,z) 


i 


/ 


dV 


f c Cl c k 
x kVA c (-r' 



(41) 


where we introduced the integration variable r = (f c /k) q. The 
Fourier components of the function are most significant along 
the directions of rapid variations of the illumination aperture. 
For a uniformly illuminated aperture, the directions are orthog¬ 
onal to the aperture edge. 

The intensity transfer function (19) is approximated by its 
value at k = 0 given by 


P«( k = 0,z) 


2 

/c 2 a 



A d 



2 

Ac(-r'). (42) 


The next non-vanishing order of (k, z) scales as k 2 and can be 
neglected, since we limit the expansions to terms linear in k. 

Substituting the transfer functions (41) and (42) in Eq. (18), 
and performing the inverse Fourier transform, we obtain the 
intensity at the detector plane (21). 
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B. Spread functions for intensity and phase 

Below we define the spread functions of phase, tilt angles and 
intensity when the detection aperture is a square. We assume 
that NA^ > NAj. In the following, to simplify the notations we 
use the dimensionless coordinates 

z / z z —y z, r / dj —y r, (43) 

where Uj = A/(2NAi), and z t - = A/(2NA 2 ) define the lateral 
and axial spatial scales, identified as the lateral resolution and 
the depth of focus of the illumination aperture, respectively 
Wavevectors are also rescaled: ka z —y k. 

Calculation of the spread functions is relatively straightfor¬ 
ward and somewhat tedious, and we provide the results without 
detailed derivations. 


The real-space representation of this function is simplified at 
the in-focus position z = 0 

g(x, 0) = j6 e~ mx ^~^ 2 ^smc (nx/2) sine ( nxf $). (46) 

The in-focus (z = 0) intensity spread functions are given by 


cos [rc(x + y)(p~ 1/2)], 

m = 1, 

cos [n(-x + y)(P — 1/2)], 

m — 2, 

cos [n(—x — y)(£ — 1/2)], 

m — 3, 

cos [ti{x — y)(P~ 1/2)], 

m = 4, 


where p(r) = (/3 2 / 2) sine (nx/ 2) sine (ny/ 2) sine (nxf$) sine {reyfi)- 
The intensity spread functions calculated at the in-focus position 
z = 0 for different quadrants are plotted in Fig. 2. 


C. The intensity spread function 

Due to the square symmetry of the illumination aperture, the 

intensity transfer functions (r) (19), where m identifies the 
quadrant of the detection lens, are given by 


g{k x ,z)g{k y ,z) + g*{—k x , z)g*{-k y , z), m=l, 

p(m) = - < s(- k x,z)g(ky,z)+g*{k x ,z)g*{-ky,z), m=2, 

4 g{-k x ,z)g{-k y ,z)+g*{k x ,z)g*(ky,z), m=3, 

g(k x ,z)g(-ky,z) +g*(-k x ,z)g*(k y ,z), m=4, 

(44) 

where where we defined a piece-wise complex-valued function 




gilnzq _ e i4:7izq 2 


e -ilKzq 2 

ilrczq 


gilnzq _ ^ 
e i4nzq{q+p) _ p 


o. 


0 < q < 1/2, 

1/2 — j6<^<0, 

-p<q<l/2-p, 

q < — j6, or q > 1/2, 

(45) 


the parameter ft = NA^/NAi > 1, and g*(q,z) is the complex 
conjugate of g(q, z). The real and imaginary parts of g(q, z) are 
plotted in Fig. 12 when NA^ = 2NAi for in-focus and out-of¬ 
focus imaging at z = Z;. 



Fig. 12. Plots of the function g(q,z) at different defocus dis¬ 
tances and in the case NA^ = 2NAi. Solid blue and red dashed 
lines depict the real and imaginary parts of g(q,z), respectively, 
calculated at z = 0; solid green and dashed magenta lines depict 
the real and imaginary parts of g(q, z) calculated at z = z z -. 


D. The spread functions for phase and light tilts 

Due to the square symmetry of the detection and illumination 

apertures, the phase transfer functions P^ (k, z) (20), where m 
identifies the quadrant of the lens, are given by 




§(k x ,z)g(k y ,z) -g*(-k x ,z)g*(-ky,z), 
g(-k x ,z)g(ky,z) - g*(k x ,z)g*(-k y ,z), 

g(-k x , z )g(-ky,z) - g*(k x ,z)g*(ky,z), 
§( k x,z)g{~ky,z) -g*(-k x ,z)g*(ky,z), 


m=l, 

m=2, 

m=3, 

m=4, 

(48) 


where the function g(q,z ) is defined in Eq. (45). 

The in-focus phase spread functions are thus given by 


sin [n(x + y)(£ — 1/2)], 

M = l, 

sin [n(—x + y)(p — 1/2)], 

m = 2, 

sin [n(—x —y)(P~ 1/2)], 

m = 3, 

sin [n(x — y)(£ — 1/2)], 

m = 4, 


(49) 


where p(r) = (j6 2 /2)sinc(7rx/2)sinc(7ry/2)sinc(7rxj6)sinc(7ryj6). 
The phase spread functions at the in-focus position z = 0 for 
different quadrants are plotted in Fig. 2. 

The in-focus tilt transfer functions are calculated as linear 
combinations of the phase spread functions according to Eq. (11) 

px _ p(4) i p( 4 ) _ p( 2 ) _ p( 3 ) 

1 Q — 1 (p i 1 (p 1 (p 1 (p / 

py = pW + pW-P^-pW. (50) 

Taking the limit z —y 0 in Eqs. (48) and substituting the result 
into (50), we obtain 

p$( k) = 2ig-(k x )g+(ky), p y e ( k) = lig- (ky)g+(k x ), (51) 

where 


g-(q) = - g(q,0)}/2, 

g+(‘l) = Il(-^O) +g{q, 0)] /2. (52) 

The functions (q) and g+(q) are the odd and even functions, 
respectively. For positive arguments, the functions are 


£-(<?) 


q, if 0 < q < 1 / 2, 

1/2, 1/2 < <7 < £-1/2, 

P~q, p — 1/2 < q < p, 

0, q > p, 


( 53 ) 
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g+(q) = * 


l-q, 

1 / 2 , 


0 < q < 1/2, 

1/2 < q <0-1/2, 
0 - 1/2 < q < 0, 


0, q > 0. 


(54) 


Similar to Eqs. (52), we introduce odd and even functions of 
the argument q within the bandwidth \q\ < 1: 


g-(q,z) 


q sine [2nzq 2 ], 

sm[27rz|(j|(l-|i)|)l 
2 nzq 


M <1/2 

1/2 < \q\ <1, 


(60) 


The functions are shown in Fig. 3 for the case ft — 2, that is when 
NA d = 2NAi. 

Real-space representations of the functions are given by 

g+(x) = Reg(x,0), g-(x) = —z lmg(x,0), (55) 


g+(q, Z ) = (l-\q\)smc[2nzq(l-\q\)}, \q\<l. (61) 

The functions g- ( q,z ) and g + (^,z) equal zero at \q\ > 1. 

The reduced phase transfer functions calculated for different 
quadrants are given by 


where Re and Im are the real and imaginary parts, respectively, 
and the function g(x, 0) is defined in Eq. (46). 

The tilt spread functions are the products of the functions 

Pe(r) =2ig-(x)g+(y), P e *(r) = 2ig-(y)g+{x). (56) 

The functions are plotted in Fig. 5 in the case /3 = 2, that is when 
NA d = 2NAi. 


E. Properties of the 3D phase spread functions 

At a defocus distance z, the phase spread functions exhibit lat¬ 
eral shearing and blurring. These effects can be traced to the 
behavior of the corresponding transfer functions g(q, z) in (45) at 
different spatial frequencies. For example, within the tilt range 
\k X/ y | < 1/2 (in dimensionless units), the transfer functions are 
proportional to the phase factor encoding lateral shearing in 

real space P^ (k, z) ~ exp [i2n(±k x d= ky)z\, where the signs 
depend on the quadrant m. 

Lateral shearing of the spread functions allows identification 
of the corresponding defocus position z. For simplicity, let us 
consider the reduced transfer functions in which one of the two 
components of the wavevector is set to zero. Depending on the 
value of the remaining second component, the transfer functions 
are related to one another. For the first and second quadrants 
we find 

Pfi\k x , 0,z) = -P^\k x ,0,z)h{k x ,z), (57) 

where 

{ e i4nzc l, \q\ < 1/2, 

e i2nzq(l-2\q\) ' 1/2 < \q\ < 0 - 1/2, (58) 

e iinz fo, 0 — 1/2 < |g| < 0. 


where 



m — 1, 

f*(k x ,z), 

m = 2, 

f*(k x ,z), 

m — 3, 

f(k x ,z), 

m = 4, 

f(ky,z), 

m — 1, 

f(k y ,z), 

m — 2, 

f*(ky,z), 

m — 3, 

f*(ky, Z ), 

m = 4, 


f(q / z)= l -^g-(q / z). 


(62) 


(63) 


(64) 



Fig. 13. Plots of Vq>(u r z) for different values of z. The lateral 
position is in units of a\ 


Certainly, (0, k y , z) = pjjP (0, k y , z). 

The reduced transfer functions of the first and fourth quad¬ 
rants are also proportional to one another 

P§\0,ky,Z) = ~P^\0,ky,Z)h{ky,Z). (59) 

We also obtain Py\k x , 0,z) = pj^\k x ,0,z). Similar relations 
can be derived for other pairs of the transfer functions. 

F. 3D spread functions when NAi = NA d 

The case of matching numerical apertures NAi = NA d is spe¬ 
cial, since, according to Eqs. (57) and (59), the pairs of reduced 
transfer functions are related to one another by the phase factors 
e i4:7izk x anc j giAnzky, respectively, for all frequencies within the 
detection bandwidth \k X/ y\ < 1 (in dimensionless units). 


The averaged spread functions (P^ (x,y,z))y 


and 


( P^ (x,y,z)) x are expressed in terms of the function V(x,z) 


<4 m Ay,z))y = < 


(Py m \x, y,z)) x = < 


V<p(x + z,z), m = 1, 
V<p(-x + z,z), m = 2, 
V<p{—x + z,z), m= 3, 
P^(x + z,z), m= 4, 

V<p(y + z,z), m = 1, 
V<p{y + z,z), m = 2, 
'Pcp(-y + z,z), m = 3, 
V 9 {-y + z,z), m= 4, 


(65) 


( 66 ) 
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where (..) Xf y is the average over x or y coordinates, and we 
defined 

V(p(u,z ) = - dq sin (2tcuc\) g-{cj,z). (67) 

The dependence of function Vq,{u,z) on the coordinate u is 
shown in Fig. 13 for different defocus distances. 

Eqs. (57) and (59) can be used to identify z. Once the defo¬ 
cus distance z is known, we can separate the effects of blurring 
and shearing. The latter is accounted for by including the cor¬ 
responding phase factors in the definitions of the tilt transfer 
functions: 

px _ e -i2nz(k x +ky) pW + e -i2nz(k x -ky) pW 
_ e i2nz(k x -k y ) p(2) _ j2nz(k x +k y ) p(?) , 
pV — e -i2nz{k x +k y ) pW + e i2nz(k x -ky) p(2) 

_ e i2nz(k x +k y ) p(3) _ e -i2nz(k x -k y ) pW_ (6g) 

The tilt transfer functions are given by 

P$(k,z) = 2ig-(k x ,z)g+(ky,z), 

Pg(Kz) = 2ig-{k y ,z)g+{k X/ z), (69) 

and can be used for phase reconstruction (38) with the resolution 
defined by the detection aperture, at least for small enough 
defocus distances. 

G. Imaging an out-of-focus phase step 

A phase step located a distance z away from the focal plane is 
characterized in frequency space by (33). 

For a system with matched illumination and detection aper¬ 
tures, NAi = NA^, the intensities detected in the first and sec¬ 
ond quadrants are obtained using the transfer functions (60): 

j(i) = ^0^ dq ^yp z ' ) cos [2nq (x + z)], 

I ® — J dq cos [27tq (x — z )\. (70) 

Fig. 10 illustrates these functions for defocus distances z = 0, 
z = 5 and z = 20 (in dimensionless units). 

H. Out-of focus imaging when NAi < NA d 

When the numerical apertures of the illumination and detection 
are not matched, then, according to Eqs. (57) and (59), the pairs 
of reduced transfer functions are related to one another by the 
phase factors e i4l7lzkx and c z4zrz ^ only within the tilt bandwidth 
\k Xf y \ < 1/2 (in dimensionless units). The phase factor relating 
the functions changes for spatial frequencies outside this range. 
In that case, phase reconstruction of out-of-focus objects remains 
possible, although with a spatial resolution reduced compared 
to the case of matched numerical apertures, at least for small 
enough defocus distance z. 

The formalism presented in Section F remains applicable pro¬ 
vided the bandwidth of the transfer functions (62), (63) and (69) 
is reduced to the tilt range \k X/ y\ < 1/2 (in dimensionless units). 
The transfer functions can be used for identification of defo¬ 
cus distance and phase reconstruction (38) with the resolution 
defined by the tilt range for relatively small defocus distances. 
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